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ABSTRACT 



The classical Levinson-Durbin linear prediction formulas for real valued input sequences 
are examined and compared to the recently proposed split- Levinson formulas. Both the 
autoregressive linear predictor model and the adaptive lattice model are used to formu- 
late the new split- Levinson algorithms. A brief introduction to the theory of symmetric 
polynomials is presented to form the basis of the new algorithms. Computer simulations 
are used to test and compare the computational accuracy of the new algorithms for AR 
filter coefficient estimation, parameter estimation for a moving average process, and 
spectral estimation of sinusoids in white noise. Research results indicate that the new 
algorithms reduce the number of real multiplications required for a k th order AR filter 
problem by one-half, and they are applicable to both the extended Prony method of 
spectral estimation and the estimation of moving average parameters. 
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I. INTRODUCTION 



A. OBJECTIVE 

The classical Levinson algorithm is known to provide solutions to real valued, linear 
systems involving Toeplitz structures. The computational cost for these solutions, is 
known to be 0(k 2 ). where k indicates the filter order. It has recently been proposed that 
the classical algorithm may be transformed into 2 simpler algorithms, using the theory 
of symmetric polynomials, and that either of these algorithms can be used to to solve for 
the predictor polynomial of order k at a reduced computational cost. [Ref. 1: p. 4701 

These new algorithms are termed the split-Levinson algorithms because their basis 
is formed from the concept of symmetric polynomials. These are not new in theory, but 
the application of the process to linear prediction is a new concept. Symmetric 
polynomials are based on the Barlett Bisection Theorem [Ref. 2 : pp. 1074-1076], where 
a system that possesses symmetry about a point, such as a Toeplitz matrix, can be de- 
composed into a symmetric and an antisymmetric part. The unique point of the theory 
is that either part may be used to solve the problem, or a combination of both pans can 
also be used in the solution. During our research we shall only consider real data se- 
quences. 

The split-Levinson case also has a lattice structure as the classical case. However, 
as will be shown, the structure of this lattice shows little resemblance to its classical 
counterpart. A derivation of a revised split lattice structure, and its recursive algorithm 
was attempted in order to represent the split lattice in a form similar to that of the 
classical structure. Although unsuccessful, the derivation procedure is presented for 
subject matter continuity. 

Computer programs have been written to implement the new algorithms and com- 
pare them to the classical algorithms. Additionally, computer programs are included to 
apply the split-Levinson algorithm for two cases, where the computational efficiency of 
the new algorithm could be of substantial benefit. These cases include the Moving Av- 
erage (MA) problem, where the parameters of a MA model must be determined from the 
given data, and an extension of the Prony method of spectral estimation, where a least 
squares estimation of the presence of sinusoids in white noise is made from the output 
data sequence. 
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This thesis compares the classical and lattice structures of the Levinson recursion 
formula given in [Ref. 3: pp. 145-167], and examines not only the formulation of the 
recursion formulas for these algorithms, but also the complexity of the computations 
and the resulting structure of each of the algorithms. 

B. THESIS ORGANIZATION 

The structure of the thesis is divided into 4 chapters, including the Introduction. In 
Chapter II we will review the classical Levinson algorithms. In the first case, the algo- 
rithm is obtained using the autocorrelation function of the input sequence, and in the 
second case, it is obtained using the forward and backward error vectors of the input 
sequence. In each case we shall establish the number of real multiplications required to 
complete a k-th order recursion of the respective algorithm. As stated, the ultimate goal 
is to establish the computational efficiency of the split- Levinson algorithm over the 
classical Levinson algorithm. Chapter III deals with the derivation of the split- Levinson 
algorithms preceded by an introductory section on symmetric polynomials. As in Chap- 
ter II, both the autocorrelation function and the lattice algorithms will be developed. 
In addition, a comparison between the computational cost of the Levinson and split- 
Levinson algorithms and an attempt to define the split lattice structure in terms similar 
to the Levinson based lattice are presented. 

In the final chapter two practical applications of the split-Levinson algorithm are 
investigated. These are: (1) the \1A parameter estimation problem, and (2) the extended 
Prony method. In case (1), the Levinson recursion used to determine a predictor coeffi- 
cient vector is replaced by the split-Levinson algorithm. A comparison between the test 
coefficients and the computed coefficients is presented. In case (2), an estimation of 
sinusoids in white noise is performed. Additionally, overall conclusions of the research 
as well as proposed topics for continued thesis research are presented. 



II. THE CLASSICAL LEVINSON ALGORITHMS 



The importance of the Lesinson algorithm in linear prediction theory is well known. 
The reason to present the algorithm in its two forms is twofold: (1) to present certain 
definitions that will be required later in the development of the split-Levinson algo- 
rithms. and (2) to detail the computational complexity of the Levinson algorithm for 
comparison purposes to the split-Levinson versions of the algorithm. In the context of 
our discussion within this thesis, \vc shall confine ourselves to the study of autoregressive 
modeling problems of real sequences as in Figure 1. (Ref. 3: p. 152] 




Figure 1. Autoregressive Model. 
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We know from linear prediction theory the augmented normal equations given by, 

R x & = r, xy (2.1) 

can optimally be solved by the Levinson algorithm, and that this algorithm can be im- 
plemented with either the autocorrelation function or the forward and backward pre- 
diction error vectors of the input sequence [Ref. 3: pp. 152-170]. 

A. THE LEVINSON ALGORITHM 

In order to examine the computational complexity of the Levinson recursion, it is 
necessary to formulate the recursive algorithm, and to determine the number of real 
multiplications and additions required to complete the algorithm. First, construct a 
Toeplitz matrix from the sequence s(t), of length N, defined as R k = 0 < /, j < A], 

where the elements of the matrix are the autocorrelation lags given by [Ref. 2: p. 646] 

A — 1 —i 

Ri= v _‘i_. X *« *(' + ') (2-2) 

i= 0 

then, the predictor vector a can be determined as a solution to the system defined by the 
matrix equations 

= of (2.3) 

where o k is the prediction error norm, and is defined as 

k 

°k = ^0 + ^ a ki Rl (2-4) 

i=\ 

It is recalled from linear prediction theory, that given a positive-definite matrix R k 
of order k + 1, the kth order a coefficient vector can be computed recursively from the 
nested Toeplitz submatrices, and their respective successive predictor vectors, a. The 
well-known Levinson recursion formula is this solution, and has the form 

a ki = a k-\,i + Pk<*k-\jc-t ' = 0,1.2 k (2.5) 

with the conditions that a k<0 = 1, and a k _ Uk = 0. The parameters, p k = a kk , are called re- 
flection coefficients, also PARCOR coefficients, because they represent the partial cor- 
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relation between the zero-th and the k-th element of the prediction vector with the effect 
of all the intermediate elements removed. [Ref. 4: p. 53] 

To construct the Levinson recursion we must use the prediction error norm re- 
lationship 



~Pk)°k - 1 



( 2 . 6 ) 



and the identity 



k - 1 

°k-\Pk~~ Rk-i a k-\,l 

i=0 



(2.7) 



to define the recursion variables. Consider the following definition as it applies to the 
Levinson recursion [Ref. 1: p. 472]. 



ft— 1 

h = Rk-i a k-\,i 



i = 0 

= °k-\Pk 



( 2 . 8 ) 



and solving for p k from Eq.(2.8), we have 

Pk = 4t W 

The error norm o k can be written in terms of the the normalizing term ). k by rewriting 
Eq. (2.6), and making a substitution from Eq. (2.9) 

= -pl)°k-\ 

= °k-\ -p k (o k _ lPk ) (2.10) 

= °k-\ ~ Pk^-k 



Combining Eqs. (2.5), (2.8), (2.9), and (2.10), we have the basis for the Levinson algo- 
rithm, and it is summarized in Table 2 of Appendix A. 

B. LEVINSON LATTICE REALIZATION 

If we are given a real sequence of signal values s(0).s(l), ..., s(N-l), and it is known 
that s(t) = 0, for — 1 > t and t > .V, then for the linear prediction problem of order k we 
find it necessary to find a set of real numbers a w , a k[ , ... , a kk that will minimize the for- 
ward and backward prediction errorvectors using a linear combination of the past signal 
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vectors. If we call the forward prediction error vector f k (t) and the backward error vector 
b k {t ), and define them in terms of the a kl coefficients. [Ref. 2: p. 646] we have 



/*(')= *(' -0 ( 2 . H ) 

(=0 

k 

(- 12 ) 

i=0 

then it turns out that the same real numbers, a kl , will provide the solution to either of 
the forward or backward prediction problems, (i.e., minimize the squared Euclidean 



norm of both f k and b k ). 

Let G k be defined as the squared norm, that is 

o* = liy 2 = f Jr* (2.13) 

From Eqs. (2.11) and (2.12) forming the first three terms of each error vector we have 
the following, 

MO) = a* 0 *(0) + a k \s( -1) + ... + a kk s{ -k) (2.14) 

/*(!) = a k0 s ( l ) + °kA°) + a k2 <~l) + ■■■■ + a k A 1 - k) (2.15) 

M 2) = a kA 2 ) + a*is(l) + a*2^(0) + ... + a»s(2 - k) (2.16) 

and, 

&*(0) = a ^(°) + a kk -A -1) + a kj<-2 s ( -2) + - + %o 5 ( -k) (2-12) 

M') = a*jts(l) + ^>-1^(0) + a kk -A -*) +••• + «*o*(l - k) (2.18) 

b k ( 2) = a k A 2 ) + + a kjt-2 s (°) + - + a kO s ( 2 ~ k) (2.19) 
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If we examine the elements of these two vectors we can see they are related in that 
each can be derived from the other by reversing the order of the a coefficients. If we form 
the Euclidean norm of each vector, |[fj| and |jbj|, we see that the k- th predictor vector 
[a w , minimizes the error norm, and |)/'*j| = ||£J. 

From [Ref. 3: pp. 156-157], we can use the Levinson algorithm to define the 
recursion formula for the forward and backward prediction errors given by 



If we let the following definition apply to the lattice version of the Levinson algo- 
rithm 



then using Eqs. (2.9), (2.10), (2.20), and (2.21) we can summarize the Levinson lattice 
algorithm in Table 3 of Appendix A. 

Even though the lattice algorithm is implemented directly from the data samples, its 
computer implementation will be more complex because of the vector manipulations 
that must occur in each iteration. The lattice structure defined by Eq. (2.20) is shown in 
Figure 2. 

In summary, we discussed the Levinson algorithm which uses the autocorrelation 
elements in its recursion and the related lattice structure which uses the input data di- 
rectly in its formulation. In terms of the computational complexity, both algorithms re- 
quire real multiplications of the order k 2 , as detailed in Tables 2 and 3 of Appendix A, 
in order to realize a k th order predictor filter. 



fk(')=fk-i(‘) + Pkbk-ii'- •). 
bk(‘) = Pkfk-M + b k-\(‘ ~ 0 



( 2 . 20 ) 



\+k-2 




( 2 . 21 ) 
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Figure 2. Lattice Prediction Error Filter. 
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III. THE SPLIT-LEVINSON ALGORITHMS 



The split- Levinson algorithms are based on the theory of symmetric and antisym- 
metric polynomials. We know that for an k.-th order real autoregressive process, the 
normal equations are 



1 

o 

’X) 

1 




1 




a 


R x R 0 Ri ... R k _ { 








0 


R 2 R\ Rq ... R k _ 2 




*2 


= 


0 


&k Rk- 1 R-k- 2 ••• ^0 




a k 




0 



or, 

Ri k - [<j, 0,0 Of 



(3.2) 



Using the Barlett Bisection Theorem [Ref. 5: pp. 1074-1076], and because of the sym- 
metry of the autocorrelation matrix, we can say that the predictor coefficient vector is 
the linear combination of a symmetric and antisymmetric predictor vector given by 

+ ( 3 - 3 ) 

The symmetric and antisymmetric vectors are defined as 



a 



(S)_ 
k ~ 



a< a) - 

*k ~ 




(3.4) 



where B represents one-half of the vector components of a*, and B' represents the re- 
versal of the vector components of B. Using Eqs. (3.3) and (3.4), we can transform the 
normal equations into 
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(3.5) 



Raf- La) 2, 0.0 c/lY 

Rtf = 0/2, 0,0 -ol2tf 

Therefore, we can see some favorable consequences of these revised normal 
equations, and their solutions. First, either the symmetric or antisymmetric form will 
give the same solution, and second, because of the symmetry of the predictor vectors, 
we need only solve for one-half of the predictor coefficients. 

Similar to the Levinson algorithm we now proceed to develop the split-Levinson 
algorithms from the input sequence autocorrelation function and the predictor error 
vectors. 

A. SPLIT-LEVINSON ALGORITHM 

The predictor polynomial a k (z) is defined as 



k 

a k( : ) = Xj a/!iZ ~‘ 

i= 0 

relative to the given Toeplitz matrix of autocorrelation lags. Denote the reverse of our 
predictor polynomial as a(z) = z"*a*(z _1 ), and the predictor polynomial has been shown 
to obey the recursion [Ref. 3: pp. 156-157] 

a k {z) = + P k z~'a(z) P- 7 ) 

and the reverse polynomial of Eq. (3.7) is 

a k (z) = z~'a k _ x (z) + p k a k _ x (z) (3.8) 

We now want to form a new polynomial from the given predictor polynomial that 
will form the basis of the split-Levinson algorithm. It is desired to show that the deter- 
mination of the coefficients of this polynomial will allow us to recover the original pre- 
dictor polynomial, and at the same time be more computationally efficient. We define 
the symmetric polynomial as P k (z ), and the antisymmetric polynomial as P[ B) {z) , and we 
desire them to be of the form [Ref. 1: p. 472] 
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(3.9) 



k 

p k (z) = y Pu z-‘ 

1 = u 

rPM-irfr - 1 

Recall from Eq. (3.4) and (3.5) that the symmetric and antisymmetric predictor coeffi- 
cients are composed of two vectors that are reverses of each other, and we will define 
these vectors so that they obey the relationships 



Pki Pkjc-i 

J*) _ _ Jfl) 
Pk,i Pk,k-i 



(3.10) 



Consider the mathematical interpretation of making the autocorrelation matrix, R k 
a singular matrix. If the reflection coefficient p k is made + l, then this corresponds to 
an element of R k making the matrix singular. For this reason we shall designate the 
symmetric and antisymmetric predictor polynomials as singular predictor polynomials 
[Ref. 2: p. 472] and from Eq. (3.9) they are defined as 

p k {z) = a*_i(z) + r'a t .i(:) 

n(a)/ \ t s -1 A / N U 

Pk W = «*-i ( z )~ z a k - iW 



Also, these singular predictor polynomials are self-reciprocal [Ref. 2: p. 472] because of 
their symmetry and may be expressed in the following forms 



P k (z) = z- k P k (z~') 
Pf = -Z~ k P k ( >"‘) 



(3.12) 



From Eq. (3.1 1) we have 



= a k _ l (z)-P£\z) 



(3.13) 



If we add Eqs. (3.7) and (3.8), and make a substitution from Eq. (3.13) we have 

M z ) + <*k( z ) = z ~ l a k _ Y {z) + Pua^iz) 4- a k _ x (z) + p k a k _ x (z) 

= (l + Pk) a k -\( z ) + 0 + Pk) z la k- i( z ) 
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where we have defined X k as 



>-k = 1 + Pk (3.15) 

In a similar fashion we can solve for the antisymmetric normalized singular predictor 
polynomial by subtracting Eqs. (3.7) and (3.8), and substituting from Eq. (3.13) we have 

a k (z)-a^ = >.fp^\z) (3.16) 

where 

4 a> =l ~Pk (3-17) 



Similar to the predictor polynomial a k (z ), we can define the singular predictor coef- 
ficient vectors for Eq. (3.9) as [Ref. 2 : p. 472] 



Pyfe “ \-PkQ'Pk\y 



n (a) — 
Pk ~ 



= O*0. Pk\’ 






(3.18) 



Since we want the split- Levinson normal equations to be of the form 



Rfa = l>*p,o o] r 


(3.19) 


R k z k = [0,0 a k ] T 


(3.20) 



then from Eq. (3.14) and (3.16) the singular predictor polynomials are 



P k (z) = 



p ( k%)= 



a k(z) , 



'• k 

*k(z) 

;(«) 

A k 



'■k 

;(«) 

A k 



(3.21) 



Since a k {z) is a polynomial formed from the predictor coefficient vector that is a solution 
to Eq. (2.3), it follows that P k (z) and Pi a) {z) are solutions to the Toeplitz system described 
by Eqs. (3.19) and (3.20). [Ref. 1: p. 472] 
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Normalizing Eqs. (3.19) and (3.20) by ). k and the split-Levinson normal equations in 
matrix form are expressed as 



If we expand the matrix expressions in Eq. (3.23), the modified error norms may be ex- 
pressed as finite sums of the predictor coefficients 



where R, is the i-th autocorrelation element of the k-th sub matrix. 

Since the symmetric and antisymmetric polynomials are closely related, we shall 
derive only the symmetric polynomial recursion equations, and then simply present the 
results for the antisymmetric case. 

The final step in the derivation is to derive a three term recursion formula for the 
symmetric polynomial. From Eq. (3.1 1) and (3.14) we have the surprising result that the 
predictor polynomial a,{z) can be obtained from a linear combination of successive sin- 
gular predictor polynomials [Ref. 2: p. 472]. First, form P k ^(z) from Eq. (3.11), and then 
eliminate a k (z) using Eq. (3.14) 



= C^.o t o, - ,T*] r 

KtpSSOP.o.o -VY 



(3.23) 



where we have defined the modified prediction error norm [Ref. 2: p. 651] 




(3.24) 



k 




im 0 

k 



(3.25) 
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(3.26) 



^*+i to = «*(*) + z"‘4(-) 

= a k (z) +z~' \_/. k P k ( z) - a*(z)] 

= (l-z- 1 K(z) + z-‘4P,(z) 

(!-z- , K(z) = ^ + iW-z-V,(z) 

If we replace k by k-1 in Eq. (3.26) above we also have 

(1 - z" Vito = P k (z) - z~ l A k _ t P fc_,(z) (3.27) 

We now form our recursion formula by mutiplying Eq. (3.1 1) by (1 - z~ l ), and use Eqs. 
(3.15), (3.26), and (3.27) to eliminate p k and all a k predictor polynomials. 

(1 - z _l )a*(z) = (1 - z _l )«*-ito + Pt(l ~ z - ')z -1 ^_i(z) 

= (1 - z _1 )«*_,(z) + p*(l - z"')[P*(z) - a*_,(z)] 

= (1 - z _1 H-ito + P k [P k (z) ~ %_,(z) - z~'P k (z) + z” Vito] (3.2S) 
= «fc-i(-)[l - r"' - p k + z~'p k ] 

= « fe _ 1 (z)C(l -z-'XI -A>fc)] 

If we now substitute for (1 - x z- x )a k (z~ x ), (1 - z- x )a k _ x (z~ x ) from Eqs. (3.27) and (3.28), and 
eliminate p k using Eq. (3.15), we can complete the derivation 

Vito - z-'A k P k (z) = [/>*(z) - z _1 Vi Vito]( 2 - ;. k ) + V. k P k (z) - P*(z)](l - z- 1 ) 

= 2 P k (z) - ?. k {z)P k (z) - 2 z~ Vi^-ito + z~'h> k -\P k -\(A 
+ >. k P k (z ) - P k (z) - z~'). k P k {z) + z~'P k (z) (3.29) 

V,W = 0 - *~V*to + z“ VA-.toP* - 2] 

= (1 + z~')P k (z) - a k z~' P k _ x (z) 

Taking the inverse Z transform of Eq. (3.29), we have the three term recursion formula 
for the singular predictor coefficients 

Pki — Pkl + Pkjc-i ~ a kPk-\Ji-i (3.30) 

where the recursion parameter a* is defined as 

«*-V,[2 - 4 ] (3-31) 

We note that x k is determined from P k {z) from Eq. (3.23), and therefore we conclude that 
the singular predictor polynomials can be recursively computed from Eq. (3.30). How- 
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ever, the recursion parameter a* is not quite in the correct form. From Fqs. (2.10), 
(3.15), and (3.31) we can alternatively compute x k from [Ref. 2: p. 473] 



T *-i 



(3.32) 



The dual relationships for the antisymmetric split- Levinson formulas can be derived 
by following a procedure similar to the one presented above. It suffices to replace the 
quantities /*, r k jx k , p h , by their antisymmetric duals, i.e., /?£>, and use the following anti- 
symmetric initial conditions. [Ref. 2: p. 649] 



/4?-o 

P \ 0 = 1 
P\\ — ~ 1 

T ( u ) _ D 

T 0 - K o 



(3.33) 



Recursive equations for the symmetric split- Levinson -algorithm are summarized in 
Table 4 of Appendix A. Examining the entries in Table 4, we see that a full iteration 
loop of the algorithm requires approximately t real multiplications. However, because 
of the symmetry of the singular predictor coefficients, we only have to perform one-half 
of these calculations. Therefore, for a k-th order filter we need to make on the order of 
k 2 /2 real multiplications. The <5 function in Table 4 is used to distinguish between even 
and odd orders of the indexing variable. 

The FORTRAN program SPLIT, in Appendix B, estimates the predictor coefficients 
using the Levinson and split-Levinson algorithms. Figure 3 is a graphical comparison 
between the known test filter coefficients of SPLIT, shown by the solid curve, and the 
filter coefficients computed by the Levinson and split-Levinson algorithms, shown by the 
dashed curve. We now undertake the derivation of the lattice form of the split-Levinson 
formulas to verify the numerical complexity of that method, and to investigate the 
symmetric and antisymmetric lattice structures compared to the Levinson lattice forms. 

B. SPLIT LATTICE ALGORITHM 

We begin the split lattice derivation by introducing the symmetric and antisymmetric 
error vectors x k , x[ a) [Ref 2: p. 648]. If we use previously established singularity concepts, 
and substitute ± 1 for p k in Eq. (2.20) for the symmetric and antisymmetric error vectors, 
x* and x?\ respectively, we have 
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Figure 3. Levinson / Split Levinson Coefficient Comparison. 

x k( { ) =fk - iW + ^-i( f “ !) ,3 3 n 

As in the split-Levinson case, we shall proceed with the derivation of the symmetric 
split lattice, and present the antisymmetric lattice results at the end of the derhation 
with any significant changes noted. 
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If we extend the singular polynomial concept to the singular predictor coefficients, 
we can start with the Levinson coefficient recursion formula and substitute ± 1 for p k 

a ki = a k-U + Pk a k-\,k-i {1.35) 



and substituting for p k 



Pki ~ a k-\,i + a k- ljc-i (3.36) 

If we write Eq. (3.34) representing the time index (t) with the subscript (i) we have an 
algorithm that is more easily adapted for computers. 

x ki=fk-u+ b k-\M-i (3-37) 



Now, comparing Eqs. (3.36) and (3.37) we have a direct correlation between the two, 
and from the split-Levinson equation for the forward error vector, we can write the dual 
split lattice equation for x K (t) 



k 

= X Pki s (t - 0 (3.38) 

/=o 

Since Eq. (3.3S) is in the form of a convolution sum we can apply Z transform theory 
to see if any inferences can be made 



X k (z) = P k (z)S(z) 



■MU 

S(z) 



- PM 



(3.39) 



From Eq. (3.39) we can conclude that the symmetric polynomial constitutes the Z 
transform of the transfer function formed from the error vector and input sequence z 
polynomials. Now we can use the previously derived split-Levinson algorithm, and in- 
verse transform it to obtain the lattice error vector recursion algorithm. 

Repeating the split-Levinson recursion we have 

^+i( z ) = (i + z "V*( z ) - v'ViM .... 

_i (3.40) 

= P k {z)+z ' P k ( z ) ~ a k z P k-\( z ) 



Multiplying Eq. (3.40) by S(z) and using Eq. (3.39) for P k {z)S{z) we have 
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S(z)P M (z) = S(z)[/>*(z) + z~'P k (z) - «*/»*_,(!)] 

= A*(r) + z-'x k (z) - «*r -l A'*_,U) 

Applying the inverse Z transform to each side of Eq. (3.41) we have the singular pre- 
dictor error vector recursion formula [Ref. 2: p. 650] 

**+](') = x k(‘) + z~'x k (t) - (3.42) 



The symmetric lattice structure given by Eq. (3.42) is shown by Figure 4. [Ref. 2: 
p. 650] 

From Eq. (3.32) we know that the recursion parameter c/. k is defined as r and 
since oc k appears in the recursion formula for the singular error vector, we need to solve 
for it. We begin with the Levinson error norm equation 



a k~ 0 Pk) a k - 1 

2a* = 2a*_ 1 (l +p*)(l - p k ) 



°k 



1 + Pk 
2a*_ l (l — Pk) — 2*1 



= 2°k-\i l -p k ) 



(3.43) 



where we have substituted r k from Eq. (3.32), and 2a k _ { (\ - p k ) is defined as ||;c A || 2 [Ref. 
2: p. 650]. 

The initial conditions for Eq. (3.42) must be examined because most cases are trivial 
except for the case of k = 0. From Eq. (3.38) we have 



*o(0 = A) 0 5 W> 



(3.44) 



and from [Ref. 2: p. 648] we define p QQ = 2. 

The antisymmetric duals are very similar to the symmetric case, and can be formed 
by replacing the symmetric variable by its antisymmetric counterpart, i.e., x\ a) (t) for 
x k (t), — p k for p k , etc. The initial conditions for the antisymmetric case are given below 



= o 

X k \t) = i(r) - s{r - 1) 



(3.45) 



The antisymmetric lattice structure given by the antisymmetric dual of Eq. (3.42) is 
shown in Figure 5. The split-Levinson lattice formulas given by Eqs. (3.37), (3.43), and 
(3.44) are summarized in Table 5 of Appendix A. If we examine the number of real 
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Figure 4. Symmetric Lattice Structure. 
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Figure 5. Antisymmetric Lattice Structure. 
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multiplications given in Table 4 and Table 5, then we can deduce the following conclu- 
sions. The split versions of the Levinson algorithm do produce a reduction in the com- 
plexity of the calculations when compared to the classical versions. The split- Levinson 
produces a reduction of one-half, and the Levinson produces a reduction of 3 2 [Ref. 2: 
p. 645]. The FORTRAN program SLATIS. Appendix C, implements the Levinson and 
split-Levinson lattice algorithms, and a graphical comparison between the known test 
coefficients, shown by the solid curve, the coefficients estimated by the Levinson lattice 
algorithm, shown by the dashed curve, and the coefficients estimated by the symmetric 
split-Levinson algorithm, shown by the dotted curve, is presented in Figure 6. 

The split lattice structures shown in Figure 4 on page 19, and in Figure 5 on page 
20 show that the classical lattice structure appears to be lost in the new algorithm. The 
distinct advantage of the original lattice structure is the modularity of the filter. In order 
to retrieve this appealing feature of the lattice filter we shall now proceed to derive a 
revised version of the split lattice structure, and see if it can have a form similar to the 
classical structure. 

C. SPLIT LATTICE REVISED STRUCTURE 

To begin, consider the second order classical lattice structure derived from 
Figure 2 on page 8. We can write the following equations for the first and second stage 
forward and backward prediction errors, 



Now substituting the equations fory|(n) and for g,(n) into the equations for f 2 ( n) and 
g 2 (n), solving for the second stage forward and backward prediction filter errors and 
taking the Z transforms yields the transfer functions 



fM) = s(ri) + p l s{n- 1) 
£,(«) = p,i(n) + i(n - 1) 

fl(n) =/i(«) + PiS\{n ~ 1) 

gi(n) = Pif\(n) + g x (n - 1) 



(3.46) 




(3.47) 



and. 




(3.48) 
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Figure 6. Levinson vs. split-Levinson Coefficient Comparison, 
which obviously yield the second order predictor transfer function A 2 (z) and its reverse 
version A 2 (z), where a 20 = l, a n = p 2 . = p, — p,p 2 = Pi(l — p 2 ). Now forming the third 

order symmetric polynomial P % { z) from our second order example, we have [Ref. 1: p. 
472] 



J’jW = A lU) + Z 3 /1 2 (z' 1 ) 

= 0 |) + (a : + a 2 ): 1 + (flj + a 2 )z 2 + agz 






(J.49) 



If we define a 0 = 1. and («, + a 2 ) = p l . then 



P i (z)= 1 +/v 1 + A- 2 + - 3 

= 1 4- p\Z 1 r : ■(/’! + 2 ' ) 

Define the following terms: 

F i(-) = 1 + P \ 2 1 
G l (z) = z'\\--p ] z) 



(3.50) 



(3.51) 



therefore, 

P } (:) = F t (z) + r 2 Gt(:) (3.52) 



Equation (3.52) defines the revised symmetric split-lattice structure, and Figure 7 gives 
a graphical representation of that structure. 




Figure 7. Proposed Split Lattice Configuration 



In order to show that the preceding classical structure can be equated with the 
symmetric polynomial derived earlier, it is now necessary to form the forward and 
backward prediction error transfer function of the new split lattice structure and com- 
pare them to the Z transform of the symmetric polynomial. Therefore, from Figure 7 
and Eq. (3.52), we can write P 3 (z) as 

r 3 (z) = I + A \z~' + z-\K , + r- 1 ) (3.53) 
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Now, if we compare Eq. (3.50) to Eq.(3.53), we see that K x = p v But, from Eq. (3.49) 
we know that p x = (a, + a 2 ), and substituting for a x and a 2 from Eq. (3.48) and (3.49) 
yields 



= P\ =Pi + PjO “Pi) (3.54) 

Let us now rederive the symmetric polynomial from the revised split lattice structure. 
From Eq. (3.50) and (3.54) we have 

Pfe) = • + (P\ - Pi - P\Pi)z~' + z" 2 [(Pi - Pi- P\Pi) + z _1 ] (3-55) 

Since we have found that the symmetric lattice can be restructured to a form similar 
to the classical lattice, the next logical step is to find a recurrence relation for the new 
lattice. Let us consider a 5 th order symmetric polynomial to determine the step-down 
recursion procedure, given by 

P s (z) = 1 + p 5 ,z _l + Pul -1 + p S2 z~ 3 + Ps^ -4 + z -5 (3.56) 



From Eq. (3.51) we have 



F 2 (z) = 1 + p n z ' + p 52 z 
G 2 (z) = z~ 2 + p 5l z~ l +p 52 



As per the observations made in Figure 7 on page 23 and Eq. (3.51), we have the lattice 
reflection coefficient for the second stage, K 2 = p sr Now reduce the order of F(z) to find 
the first stage reflection coefficient using the standard inverse Levinson recursion [Ref. 
4: pp. 156-157]. 



FM = 



F 2 (z) - K 2 G 2 (z) 

\-k 2 2 



= 1 + 



Ps\ 

I + A 2 



(3.58) 



therefore 



Ai = 



Psi 

l + K 2 



(3.59) 



Now rewriting the equations for F x (z) and F 2 (z) using the derived reflection coefficients, 
we have 
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(3.60) 



F, (z) = 1 + A',z 1 

G, (s) + A, 

F 2 (z) = 1 + A|( 1 + A 2 )z * + A 2 z 
G 2 (z) = z - 2 + A',(I + A 2 )z - 1 + A' 2 



and we know that 



A s (z) = F 2 (z) + z- 3 G 2 (z) 

= 1 + A', ( 1 + A 2 )z _1 + K 2 z~ 2 + A' 2 z -3 (3.61) 

+ A j ( I + A 2 )z + z 

Equating the terms in Eq. (3.56) and (3.61), we have p sl — K x {\ K 2 ) and p S2 = K 2 . 
Knowing the values of K x and K 2 , we now can form a two stage symmetric lattice similar 
to that in Figure 7 on page 23 to implement P s {z) . However, we are interested to find 
if we can recursively obtain the lower orders P 4 (z), P 3 (z), etc. or the higher orders 
P 6 {z), /\(z), etc. from P s (z). In an attempt to form P 6 {z), we use the standard forward 
recursion [Ref. 3: pp. 156-157] to obtain 

F 3 (z) = F 2 (z) + K 3 z-'G 2 (z) 

= 1 + A,(l + A 2 )z _1 + K 2 z~ 2 + K 2 K 3 z~' (3.62) 

+ A,A 3 (1 ■+■ A 2 )z 4- A 3 z 



and 



G 3 (z) = z“ 3 + A,(l + K 2 )z~ 2 

+ K 2 z~' + A 2 A 3 z -2 (3.63) 

+ A,A 3 (1 + K 2 )z~' + A 3 

Forming the symmetric polynomial P 6 {z) from Eqs. (3.62) and (3.63) we have 

P 6 (z) = F 3 (z) + z _ 3 G 3 (z) 

= 1 + (A', + A,A 2 + A 2 A 3 )z _1 

+ (A 2 + K 3 K 3 + A, A 2 A' 3 )z -2 + K 3 z~ 3 (3.64) 

+ A 3 z -3 + (A 2 + A, A 3 + A,A 2 A 3 )z - 4 
+ (Aj + A, A 2 + A 2 A 3 )z + z 

Comparing Eq. (3.64) to the symmetric form of P 6 (z), we have 
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Pb\ - (^1 + ^ 1^2 + ^2^3) 



P5\ P5\P52 

1 + /? 52 1 + Ps2 



+ —PS2P63 



(3.65) 



The one-half enters into Eq. (3.65) because we know that for even orders the polynomial 



coefficients are symmetric about the center element, and they must be shared in the 
matrix equations. We shall now expand Eq. (3.65) to attempt to develop a recursive 
algorithm for the sixth order coefficients from the fifth order coefficients. Expanding Eq. 
(3.65) we have 



wfiere the term in brackets is an expansion of the coefficient recursion formula, Eq. 
(3.30), for p 6 3 . Collecting terms we have 



From Eq. (3.68) we can observe that the a k recursion parameter and the number 
of previous coefficients required to be known are increasing in order, and it appears that 
a simple recursive algorithm based on the above approach is not possible. Note that 
although the new lattice structure does not appear to be order recursive, we can express 
a given order symmetric or antisymmetric lattice structure in a more conventional form. 

In summary, we know that the split- Levinson algorithm is a viable replacement for 
the classical algorithm because of its computational efficiency. We have studied both 
autocorrelation and data (or lattice) realizations of the split-Levinson algorithm. An at- 
tempt to derive a recursive split lattice algorithm yielded a classical-like lattice structure, 
but it is not recursive in order. Further investigation is necessary in this direction. We 
now need to test the split-Levinson algorithms on some signal processing applications. 




(1 + p5l)P6\ ~ ^SlU + Psi) + 0 + Psi) 2 P52LP52 + Ps 2 - a 5^423 

oc s 

Pb\ ~Ps\ + P52IP52 2 ~ 



(3.67) 



Substituting for p 52 from Eq. (3.30) 



P52 ~ P*2 + PAl ~ (X 4P3\ 



(3.68) 
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IV. APPLICATIONS OF THE SPLIT-LEVINSON ALGORITHM 



In this chapter, we apply the split-Levinson algorithm in (1) the MA parameter es- 
timation problem, and (2) the extended Prony method of spectral line estimation. Before 
we take up these two applications we examine the algorithm's usefulness if the original 
filter has coefficient symmetry, i.e., the impulse response of a linear phase FIR filter. 

A. HANKEL AND TOEPLITZ MATRICES 

In previous derivations we have assumed the FIR filter equation to be non- 
symmetric. Let us now investigate the problem where the filter equation is symmetric, 
i.e., of the form 



y(n) = -f a { s(n — 1) 4- a 2 s(n — 2) + -- 

+ tf /c _ 1 s(l) + M"-*) 

By definition, a symmetric polynomial is self-reciprocal, that is 



(4.1) 



(4.2) 



Therefore, from the Levinson algorithm, predictor polynomials are known to obey the 
recurrence relation 



a k {z) = a k _\{z) + p k z X a k _ x (z) (4.3) 

and in our special case we have 

a k (z) = a*_,(z) + P*z _1 a*_,(z) = a(z) 

-(l + p*r“ l )a*_i(r) 

In order to formulate a set of equations similar to the split-Levinson, it is necessary 
to derive the normal equations for our special case, and compare them to the standard 
equations, in order to develop the recursive algorithms. Since the predictor coefficients 
in a recursive algorithm produce estimates of s(n) as the algorithm steps through its re- 
cursive steps, denote this estimate as s(n). In vector form, we then have 
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a \ 

a 2 



s(/i|«-l)~ - [s(n- l)s(*-2)...s(0)] 






a 2 

a \ 

l 



(4.5) 



To derive the normal equations we must find the minimum mean squared error (VISE) 
from the equation for the error, 

e(n) — s(«) — J(«|/f - 1) (4.6) 

The minimum mean squared error is found by squaring the error term, and then differ- 
entiating the squared term with respect to the given a k vector. Combining these two ev- 
olutions we have the following equations, 

MSE = J 

= £|> 2 («)] (4.7) 

= £[(*(«) -(J(„)|„-l)) 2 ] 

To obtain the normal equations, we carry out the following steps: 

77 = 0 = -2£[j(«)s [_,] -2£[sJ_,s Jt _ 1 ] 

£[^)s A _,s[_,]a -£[</»)**_,] ( 4 - 8 ) 

/{(*-'> a = rf-') 

From the split-Levinson recursion formulas we know that the singular symmetric 
polynomial, P k {z ), is defmed as the following, 

p &) + ;"V,M (4.9) 

Since our predictor polynomial is symmetric, it is a reasonable question to ask if sym- 
metric polynomials also obey this recursive relationship. It is noted as an immediate 
consequence of the recursive problem, all of the preceding polynomials will also be 
symmetric. Therefore, we check the singular predictor polynomial recursion to see if it 
holds when the original polynomial is itself symmetric 
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a k {z) = a k ~ x (z) + z (z) 

= ^_l(z) + z *Z *«/ c _ 1 (z) 

= l 4- ajZ -1 + a 2 z~ 2 4- ... + z~ (/c_1) + z -/c [l 4- a t z _1 4- a 2 z~ 2 4- ... 4- z~ k ~\ 
= l 4- (1 4- ci\ )z 4- ( ci\ 4- <3|)z 4- ... 4- z 



(4.10) 



Now that we see that the Levinson recursive equation holds for a symmetric 
polynomial, we derive the recursive relationships for our polynomial using what is 
known from the split- Levinson equations. We have defined the symmetric polynomial 
P k {z) to be a normalized combination of a non symmetric polynomial, a k (z), and its re- 
ciprocal image, a k (z) in the form of, 

/ k P k {z) = a k {z) + a k (z) (4.11) 



By direct substitution it is a trivial matter to show that this relationship also holds for 
a symmetric polynomial, a k {z). 

In order to develop the recursion for the symmetric polynomial, it is necessary 7 to 
express the desired linear predictor in terms of the previous two predictors. To this end 
use Eq. (4.10), to form a*. x (z), and substitute from Eq. (4.1 1) to perform this task. 

a k (z) = a k _ x (z ) + z _1 4-iU') 

a *+i(‘) = a kU) + *"'«*(*) 12 

= a k (z) + z-'V. kak (z)-a k (z)-] 

= «*(-’)( 1 - z _l ) + z _l 



Solving for (1 — z~ l )a k (z), and forming the quantity (1 — z _1 )a A _j(z) we have 

(1 - z-’H(z) = a k+x (z) - z-'l k a k (z) 

_i _i (4.13) 

(1 - z )a k _ x {z) = a k {z) - z ). k _ x a k _ x {z ) 

These relationships will now allow us to form the three term recursion for the given 
symmetric polynomial from Eqs. (4.3), (4.11), (4.12), and (4.13) 

a k {z) = a*_,(z) + P k a k _i(z) [4) 

a M (z) = (1 + z~')a k (z) - a k z~'a k _ x (z) 



where we have defined a k 






(4.15) 
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From Eq. (4.14) we can see that the coefFicient recursion formula is the same form 
as Eq. (3.29), and we can deduce that the split- Levinson algorithms will work equally 
well for symmetric polynomials that describe unknown filters as it does for polynomials 
that are not symmetric. 

B. FIR MOVING AVERAGE PARAMETER ESTIMATION 

If we consider an FIR filter with an input sequence given by 
s r = [i(«)i(« — l)...i(rt - m)], and an output y(n) given by 

M 

yi n ) = ~ ') ( 4 - 16 ) 

i=0 

then we can develop the necessary equations to estimate the moving average parameters, 
and solve for the FIR filter coefficients. The algorithm to estimate the predictor coeffi- 
cients can be defined as follows: 

Let the three predictions, }y(n), sj(ri), and 5"(/i), represent the m-th order predictions 
of the forward estimate of y, the forward estimate of s, and the backward estimate of s 
respectively. [Ref. 6: p. 1 ] 



A/n 

yf 



m 

(n) = s(n - i) 

t = 0 



(4.17) 



m 

*/(«) = Yj C ‘ " 0 



1=0 



(4.18) 



m 

s”(n — m) = y^- s(n — m + i ), 

i= o 



(4.19) 



where the coefficient vectors are defined as, 
b r = [1 — b 0 

c r = [0 1 -c,...-c m ] (4.20) 

d r =[0 -c m -c„_, ... -cj 

Without going through all the details for the MA parmeter recursions, we can 
summarize the recursion formulas for the predictor coefficients as [Ref. 6: p. 2] 
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(4.21) 



Notice that the predictor vector d" is not included in the preceding definitions since it is 
the reverse c m 

If we examine Eq. (4.21) we see that the recursive relationship for c m is a statement 
of the Levinson recursion, since A.'?, and K™ are the m-th order reflection coefficients. 
Therefore, we can apply the split- Levinson algorithm to solve for c m , form d" , and 
recursively determine b w . Finally, from the theory of Moving Average processes, b m = a m 

The FORTRAN program MAV1, in Appendix D, uses a 25-th order FIR test filter 
to to generate a test sequence, and the results are given in Table 6 of Appendix A. 

Figure S is a graphical comparison between the known test coefficients, shown by 
the solid curve, and the computed filter coefficients, shown by the dashed curve. The 
curves are fitted to the linear magnitudes of the coefficients by interpolating spline 
techniques. 

C. EXTENDED PRONY METHOD 

The estimation of the existence of sinusoids in the presence of noise is a common 
occurrence in signal processing applications. In this simulation, we will show that the 
split-Levinson algorithm can be directly implemented in the process to determine the 
approximate frequencies. The noise is zero mean, unit variance, and white in nature. 

In this application of the split-Levinson algorithm, we attempt to approximately fit 
p exponentials to a data set of N samples. We have the constraint that .Y > 2 p, and a 
least squares estimation procedure is used. We begin by defining the estimate of our set 
of data samples. [Ref. 7: p. 1404] 



where b m — A m exp(/0J/2, and z m = exp(/27r/ m Ar). The z m 's are roots of unit modulus with 
arbitrary frequency and occur in complex conjugate pairs as long as f m 0 or l/2Af 
Therefore, in order to solve for the p sinusoids, we must solve for the roots of the fol- 
lowing polynomial. [Ref. 7: p. 1406] 



p 




(4.22) 
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Figure 8. MA Coefficient Comparison. 

= y a k z2p ~ k ~ 0 (4.23) 

k=Q 



The roots can be of unit modulus, and occur in complex conjugate conjugate pairs if 
we constrain the polynomial coefficients to be symmetric about the center element of the 
polynomial [Ref. 7: p. 1407]. This is an exact ocurrence in the symmetric and anti sym- 
metric polynomials of the split-Levinson algorithm, as long as the order of the 
polynomial is even, that is 



a 



T 



[■ 



a \ 




(4.24) 
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